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We study equal and unequal-mass neutron star mergers by means of new numerical relativity sim¬ 
ulations in which the general relativistic hydrodynamics solver employs an algorithm that guarantees 
mass conservation across the rehnement levels of the computational mesh. We consider eight binary 
configurations with total mass M = 2.7 Mq, mass-ratios g = 1 and q = 1.16, and four different 
equations of state (EOSs), and one configuration with a stiff EOS, M = 2.5Mq and q — 1.5, which 
is the largest mass ratio simulated in numerical relativity to date. We focus on the post-merger 
dynamics and study the merger remnant, dynamical ejecta and the postmerger gravitational wave 
spectrum. Although most of the merger remnant are a hypermassive neutron star collapsing to a 
black hole-l-disk system on dynamical timescales, stiff EOSs can eventually produce a stable massive 
neutron star. During the merger process and on very short timescales, about ~ 10“® — 10~^ Mq of 
material become unbound with kinetic energies ~ 10®°erg. Ejecta are mostly emitted around the 
orbital plane; and favored by large mass ratios and softer EOS. The postmerger wave spectrum is 
mainly characterized by the non-axisymmetric oscillations of the remnant neutron star. The stiff 
EOS configuration consisting of a I. 5 M 0 and a I.OMq neutron star, simulated here for the first 
time, shows a rather peculiar dynamics. During merger the companion star is very deformed; about 
~ O.O 3 M 0 of rest-mass becomes unbound from the tidal tail due to the torque generated by the two- 
core inner structure. The merger remnant is a stable neutron star surrounded by a massive accretion 
disk of rest-mass ~ O.SMq. This and similar configurations might be particularly interesting for 
electromagnetic counterparts. Comparing results obtained with and without the conservative mesh 
refinement algorithm, we End that post-merger simulations can be aSected by systematic errors if 
mass conservation is not enforced in the mesh refinement strategy. However, mass conservation also 
depends on grid details and on the artificial atmosphere setup; the latter are particularly significant 
in the computation of the dynamical ejecta. 

PACS numbers: 04.25.D-, 04.30.Db, 95.30.Sf, 95.30.Lz, 97.60.Jd 98.62.Mw 


I. INTRODUCTION 

Binary neutron star (BNS) mergers are extreme events 
associated to a variety of observable phenomena in the 
gravitational and electromagnetic spectra, e.g. m- 
BNS coalescence is primarily driven by the emission of 
gravitational waves (GWs). Indirect evidence for GWs 
has been indeed inferred by radio observation of double 
pulsars msi, but a direct detection of GWs is still pend¬ 
ing. The GW signal emitted during the last minutes of 
the coalescence and merger is in the band of ground- 
based laser interferometer network made of LIGO [5] 
and Virgo m- Within the next years, this network will 
start to operate at sensitivities where ^ 0.4 — 400 detec¬ 
tions per year are expected mm- Several electromag¬ 
netic counterparts are expected both during and follow¬ 
ing BNS mergers; joint observations of the gravitational 
and electromagnetic emissions will maximize the scien¬ 
tific returns m- Neutron star mergers are usually asso¬ 
ciated to short-gamma ray burst (and afterglows) [3 [T3] . 
Although the precise injection mechanism has not been 
clearly identified, BNSs remain the most plausible trig¬ 
gers of these powerful emissions. Dynamical ejecta from 


BNS are currently the most plausible site of origin of 
heavy nuclei (A > 140) rapid neutron-capture process 
[HHIII- The radioactive decay of some of these newly 
produced heavy elements is likely to lead to strong elec¬ 
tromagnetic transients called kilonova (or macronova) 
events [T51120) . Finally, a large amount of energy is re¬ 
leased in neutrinos, produced by the merger remnant ei¬ 
ther via shocks EBins and neutron-rich outflows |23) . 
or, at lower energies, in the hot dense regions of the hy¬ 
permassive neutron star (HMNS) [231 US]. However, the 
steep energy dependence of neutrinos of the interaction 
cross sections and their moderate energies (^ 20 MeV) 
make them hard to detect. 

Modeling BNS mergers requires relativistic hydrody¬ 
namics simulations in dynamical spacetimes, i.e. the 
solution of the full set of Einstein’s field equations. 
General relativistic BNS simulations are typically per¬ 
formed in the framework of 3-1-1 numerical relativity us¬ 
ing Cartesian-grids, finite volume methods, and explicit 
time evolutions, see |26j for a review. A crucial ingredient 
in such numerical setups is the use of adaptive mesh re¬ 
finement (AMR), in particular the methods of [27], which 
were implemented for various applications in numerical 
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FIG. 1: Conservative mesh refinement for the ID advection 
equation. The plot compares the mass conservation for a 
discontinuous profile flowing into two refinement levels in the 
two cases in which the mass correction is applied or not. Top: 
mass evolution on the three levels for the non-conservative 
(dashed lines) and the conservative method (solid). Bottom: 
evolution of mass relative error on the the coarsest level for the 
non-conservative (dashed lines) and the conservative method 
(solid). The solid red line is at round-off level. 


relativity [28H32]. Nested Cartesian boxes with 2:1 grid 
spacing refinement and Berger-Oliger time stepping [27] 
proved to be a robust and stable solution for the compu¬ 
tation of black hole and neutron star mergers [SSl- 

[39] as well as rotational collapse of neutron stars [40ti42] 
or massive stars gam]. 

One of the main problems in the simulation of hydro- 
dynamical flows with finite volume methods and AMR 
techniques is to preserve global conservation of mass and 
other quantities, especially in the presence of shocks, 
contact discontinuities and large gradients. In a semi¬ 
nal paper, Berger and Colella have proposed a refluxing 
scheme which guarantees conservation across refinement 
level gSj- Essentially, the refluxing scheme enforces the 
fluxes in and out across a coarse/fine cell boundary to be 
the same to round-off level. Algorithmically it consists in 
a correction step applied to the solution at certain grid 
points after each time step. 

The importance of a conservative mesh refinement is 
illustrated in Fig. for the simplest case of the ID ad¬ 
vection equation, {dt+dx)u{t, x) = 0 with x € [—4,4] and 
discontinuous initial data, u{0,x) = 1 for a: G [—2,—0.2], 
u(0,a:) = 0 otherwise. We employ a grid composed of 
3 fixed levels I = 0,1,2 with n = 800 grid points, cen¬ 
tered around a: = 0. The coarse level has grid spacing 
Hq = 0.01 and the others are successively refined by fac¬ 
tors two, hi+i/hi = 1/2. Time evolution is performed 


with a 4th order Runge-Kutta and the Berger-Oliger 
method; fluxes are computed with a linear reconstruc¬ 
tion using the Van Leer MC2 limiter. Figure shows the 
evolution of the mass of the system on each refinement 
level, which without mesh refine¬ 

ment (unigrid) is conserved to round off precision. Using 
mesh refinement, one observes that every time the mass 
flows in/out a refinement level (see e.g. t ~ 1,1.2, 2.2,..., 
top panel), a mass violation takes place (bottom panel). 
Notably, the mass either decreases or increases in a way 
that depends on the scheme’s truncation error and that 
in general is not predictable. Also, after the wave has left 
all inner refinement levels an error in the mass on 1 = 0 is 
still present. Instead, using the Berger and Colella cor¬ 
rection step, mass conservation is verified at the round-off 
error, exactly as in the unigrid case. 

Conservative AMR schemes have been introduced in 
numerical relativity only very recently im 05]. They 
have been used to simulate eccentric mergers of both 
black holes - neutron star and double neutron star sys¬ 
tems, including the merger remnant, post-merger disks 
and ejecta 07105]. Also, they have been employed in 
massive star and core-collapse supernovae evolutions in 
general relativity [49li5T] . Recent studies of rotating neu¬ 
tron star collapse to black hole greatly benefit of the use 
of conservative AMR, and allowed an accurate calcula¬ 
tion of the gravitational wave signal [U 02] and a lo¬ 
cal comparison of the end state with black hole space- 
times 02]. 

In binary simulations one expects that conservative 
AMR can significantly improve numerical relativity simu¬ 
lations, especially simulations of the merger remnant. A 
direct comparison of the performance of a conservative 
mesh refinement algorithm in coalescing BNS systems is 
presently missing. In the context of spinning equal-mass 
quasi circular mergers, we have pointed out that the sim¬ 
ulation of the hypermassive neutron star is sensitive to 
the mesh boxes size and their extension [53]. The latter 
factors influence mass conservation (for a fixed resolu¬ 
tion), and a conservative scheme is desirable. Another 
potentially important application of conservative AMR 
is the simulations of low-density material in postmerger 
accretion disks and dynamical ejecta. Ejecta have densi¬ 
ties several orders of magnitude smaller than the typical 
neutron star maximum densities; thus, their calculation 
employing grid-based codes is very challenging. Dynam¬ 
ical ejecta in full general relativistic BNS merger simu¬ 
lations have been previously studied only in jM] |SS] in 
more detail. Those works do not employ a conservative 
AMR strategy, thus the accuracy of the result can be, in 
principle, seriously compromised. 

The purpose of this paper is threefold. 

First, we present our implementation of a conservative 
AMR algorithm and present a set of single star spacetime 
evolutions in which we assess the performances of the al¬ 
gorithm. We focus on the evolution of different single 
star spacetimes since such tests (i) received little atten¬ 
tion in the literature; (ii) are computationally relatively 
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cheap; (iii) are highly nontrivial and preparatory cases 
for the application of the code to BNs evolutions. 

Second, we apply our upgraded code to the study 
of equal-mass (mass ratio q = 1 ) and unequal-mass 
(g = 1.16) BNS system described by various equations 
of state (EOS). We directly compare results obtained 
with and without the conservative AMR. We focus on 
the postmerger dynamics and investigate the physical 
properties of the remnant. In particular we study as a 
function of the EOS and the mass ratio the following 
properties: (i) the merger outcome; (ii) mass and kinetic 
energy of the dynamical ejecta; (iii) GW spectra. 

Third, we consider for the first time the evolution of 
a BNS system with a stiff EOS and mass ratio q = 1.5 
(total mass M = 2.5Mq). This binary has the largest 
mass ratio simulated so far (see also [56]). The partic¬ 
ular combination of EOS, g, and total mass properties 
lead to a peculiar merger remnant composed of a stable 
massive neutron star surrounded by an extended, mas¬ 
sive accretion disk. Also, the binary configuration fa¬ 
vors mass ejection during merger. These kind of binary 
configurations are possible and might be particular rel¬ 
evant for electromagnetic counterparts. However, they 
have received little attention in numerical relativity, al¬ 
though some recent observations are in favor for a stiff 

EOS |571|SH|. 


The article is structured as follows. After a brief review 
of the equations (Sec. 0 , we present our numerical strat¬ 
egy in Sec. |III| focusing on the novel implementation of 
the conservative mesh refinement. Section describes 
the main quantities employed for the analysis of our BNS 
simulations. In Sec.jVjwe describe a variety of single star 
tests in which the performance of the conservative AMR 
is investigated for different combinations of the relevant 
parameters of the simulations (restriction and prolon¬ 
gation operators and artificial atmosphere parameters.) 
Section VI summarizes the BNS conhgurations and the 
grid setup used for evolutions. In Sec. |VII| we apply the 
new algorithm and evolve 16 BNS systems with mass ra¬ 
tios g = 1 and g = 1.16 and different EOS. In Sec. |VIII| 
we consider a BNS with g = 1.5, total mass M = 2.5Mq, 
and the stiff equation of state MSlb. Einally, the conclu¬ 
sions are presented in Sec. m Throughout this article, 
geometrical units c = G = Mq = 1 are employed unless 
otherwise stated. At some places units of Mq are given 
explicitly for clarity. 


II. SUMMARY OF THE EQUATIONS 

Let us summarize briefly the most important equations 
employed in this work. We work with the 3-1-1 formalism 
(e.g. jSSj) and indicate with 7 ^ the 3-metric, and with a 
and the lapse and shift vector. 

General-relativistic hydrodynamics (GRHD) equations 
are solved in conservative form. 


( 1 ) 


TABLE I: Piecewise polytropic EOS parameters. For all 
EOSs we use a crust with Ko = 8.94746 • 10“^ and Fq = 
1.35692, and pi = 8.11940 • 10“'‘;p2 = 1.62003 • 10"T 
Columns: EOS, the density were the crust ends, the poly¬ 
tropic exponents for the individual pieces Fi, the maximum 
supported gravitational mass Mmax, the maximum supported 
baryonic mass, and the maximum adiabatic speed of sound 
Cs max within the maximum stable neutron star configuration. 


EOS 

po • 10-^ 

Fi 

Fa 

Fa 

^^max 

max 

Cs max 

MSlb 

1.84128 

3.456 

3.011 

1.425 

2.76 

3.35 

0.99 

MSI 

1.52560 

3.224 3.033 

1.325 

2.77 

3.35 

1.00 

H4 

1.43830 

2.909 

2.246 

2.144 

2.03 

2.33 

0.72 

ALF2 

3.15535 

4.070 

2.411 

1.890 

1.99 

2.32 

0.65 

SLy 

2.36900 

3.005 

2.988 

2.851 

2.06 

2.46 

1.00 


with q = Si,T) being the vector of the Eulerian 

conservative variables defined in terms of the primitive 
variables as, 

D = Wp, = W^phvi, T = (W^ph -p)-D. (2) 

The primitive variables are the rest-mass density p, the 
pressure p, the specihc internal energy e, and the 3- 
velocity uL Additionally, we dehne the Lorentz factor 
W = l/y/l — ViV\ the enthalpy h = 1 + € + p/p, and the 
determinant of the 3-metric 7 . On the right-hand-side of 
Eq. 0 one has the divergence of the fluxes and source 
terms depending on the metric, metrics first derivatives 
and fluid variables. We stress that only the first equation 
of 0 is a “strict” conservation law, 

= 0 , (3) 


in the sense that the source term is zero and a conserved 
quantity can be associated: the rest-mass Mf,. We refer 
to [Hiniin] for an extensive discussion of these equations. 

The PDE system in 0 is closed by an equation of 
state (EOS) in the form p = P{p,e). A simple EOS 
is the T-law P{p, e) = (T — l)pe, or its barotropic ver¬ 
sion P{p) = Kp^ (polytropic EOS). Several barotropic 
- zero-temperature EOS developed to describe neutron 
star matter can be fit with piecewise polytropic models, 
and efficiently used in simulations. In our work we em¬ 
ploy four segment fitting models following the construc¬ 
tion of |62j . Each segment is given by a certain rest-mass 
density interval pi < p < pi+i; the pressure is then cal¬ 
culated as P{p) = Kip^' where the polytropic constants 
Ki are determined by demanding continuity of P{p) at 
the interfaces, Kip^' = The parameters of 

our EOS are reported in Tab. [0 notice that we specify 
Po in our units. Thermal effects are simulated with an 
additive thermal contribution in the pressure in a T-law 
form, Pth = {Tth-^)pe, with Tth = 1-75, see [371IM1IS3]- 

The Einstein equations are written in 3-1-1 form, ei¬ 
ther as the BSSN [64H66] or the Z4c |67| |M| system. In 
the gauge sector, we use the H-log-slicing condition [51] 
for the lapse and the Gamma driver shift [ 711111 ]. The 


dtq = -dif + s , 
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fundamental role of this gauge in the numerical simu¬ 
lation of the gravitational collapse and singularity for¬ 
mation/evolution was investigated in different physical 
scenarios [1117M71]- 


III. NUMERICAL METHOD 


In this work we use the numerical relativity methods 
implemented in the BAM code [32J [331 ISHl ES]- Our 
new implementation of the conservative mesh refinement 
for hydrodynamics fields is based on the Berger-Colella 
method [iS] and follows [33]; we describe it in detail in 
this section. 


A. Computational Grid 


The computational grid is made of a hierarchy of cell- 
centered nested Cartesian grids. The hierarchy consists 
of L levels of rehnement labeled by Z = 0,..., L — 1. A re¬ 
finement level I has one or more Cartesian grids with con¬ 
stant grid spacing hi and n points per direction. The grid 
spacing in each refinement level is refined according to 
hi = hol2h The grids are properly nested in such a way 
that the coordinate extent of any grid at level /, Z > 0, is 
completely covered by the grids at level Z — 1. Some of the 
mesh refinement levels Z > Z""' can be dynamically moved 
and adapted during the time evolution according to the 
technique of “moving boxes”, e.g. [311 [351 ES]- BAM’s 
grid can be further extended in the wave zone using a 
multipatch “cubed-sphere” as described in [68l ElEa]. 
Every refinement level has buffer zones populated by in¬ 
terpolation; interpolation from the parent (coarse) to the 
child (fine) level is the prolongation (P) operation, the 
one from the fine to the coarse level is the restriction (R) 
operation. For metric variables these operations are per¬ 
formed with sixth order Lagrangian operators. Spatial 
interpolation of matter variables is discussed below. 

The grid variables are evolved in time with the method 
of lines, using an explicit fourth order Runge-Kutta and 
employing the Berger-Oliger (BO) algorithm [37]. For 
efficiency, we typically use only six buffer zones and per¬ 
form a linear interpolation in time to update the buffer 
zones during the Runge-Kutta step, see (33] for more 
details. A Courant-Friedrich-Lewy factor of 0.25 is em¬ 
ployed in all runs, if not stated differently. Standard 
finite differencing 4th order stencils are employed for the 
spatial derivatives of the metric. GRHD is solved by 
means of a high-resolution-shock-capturing method [33] 
based on primitive reconstruction and the Local-Lax- 
Friedrich’s (LLF) central scheme for the numerical fluxes. 
Primitive reconstruction is performed with the 5th order 
WENO scheme of jBD] as in |5T] . 


B. Conservative mesh refinement 


Let us review the main idea of the new conservative 
AMR algorithm implemented in the BAM code. With¬ 
out loss of generality we restrict the presentation to the 
first equation of Q, and to the flux in the a:-direction 
only. Although the algorithm is applied to all the fluid 
variables, Eq. the U-equation, is the only one which 
is a strict conservation law. Directions different from the 
X direction are treated in a similar way. 

The discrete model equation reads. 


T^n+l _ p.n _ 


At 

Ax 


{fI 


+ 1/2,j,k F^_i^2,j,k^ ( 4 ) 


where j k denotes the x-component of the numeri¬ 

cal flux across the cell face {i + l/2,j,k) (boundary of cell 
{i,j, k) and {i+l,j, fc)), Ax = h,n denotes the time level, 
and At the time step. Consider the model equation on 
two sequential levels of refinement with hi+i/hi — 1/2, 
and on cells at the boundary of refinement Z -|-1. Mass vi¬ 
olation happens during a BO step, because: (i) the buffer 
zones of level Z -I- 1 are set by prolongation (P) from level 
Z; (ii) the prolongation carries a certain truncation error, 
so the fluxes on Z -I- 1 differ from those on Z; (iii) after re¬ 
striction (R) from level Z -|-1, the solution on level Z is not 
consistent with the fluxes on Z. The process is illustrated 
in Fig. 

After the time step At, the changes of the vari¬ 

able on level Z due to the flux going through the cell 
face (q + l/2,ji,ki) is given by 




il + l/2,ji,ki 


(t). 


(5) 


After level Z, level Z -I- 1 advances by two At/2 time steps 
and one has 


6D 


h+i) 

{i,j,k) 


{t At) 


At/2 

Ax/2 

At/2 

Ax/2 




( 6 ) 




(t -h At/2) . 


In general, these two changes are different at truncation 
error level. Similarly the mass flows across the face are 
different, ^ SMh\ and, after restriction, the 

mass conservation is violated in a way oc . 

The original Berger-Colella algorithm corrects the so¬ 
lution at level Z after the refinement level Z -|- 1 has 
completed its time step and both levels are time- 
aligned [33]. The correction (C) operation is i—>■ 

Dh) -(- At/Ax 6Fh\ where is a flux correction 

stored on the cell face. First, SF^l is initialized with 
—F /_/^^2 ji ki before advancing in time level Z -|- 1. Then, 
during each time step of level Z -|-1, it receives and sums 

ix) 

up the contributions F . , (two contributions 

in our example). The C step guarantees consistency of 
the fluxes. East et al. [35] proposed to store the mass 
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FIG. 2: Sketch of the mesh-refinement. We focus on the buffer region along the positive a;-direction. Light red cells refer to 
the buffer region between level I and level I + 1. We employ six buffer points in level I + 1. Prolongation (P) and correction (C) 
steps take place in this region. The parent cell is visualized by the blue bounding box, while the child cells are colored dark 
red. The firrxes across the physical domain and the refinement buffer zone are visualized with arrows. The parent cell (level 1) 
receives the correction after level Z -I- 1 has been evolved. 


correction rather than 5F^^\ and perform the cor¬ 
rection as I—> jV^ '^ where is the cell 

volume |46j . This method is simpler and has the advan¬ 
tage of using grid variables defined on cell centers instead 
of faces. We follow this approach. 

Our implementation is as follows: 


1. We introduce a mask to label the cells involved in 
the C step. These are the innermost buffer points 
of level I + 1 (red in Fig. and the correspond¬ 
ing parent cells (blue in Fig. [^. The mask also 
stores the information about the box face, i.e. one 
of the possibilities {±x,±y,±z). The mask has to 
be recomputed after each regridding step. 

2. After each evolution step we store the mass change 
of the parent cells 

, (7) 

and, similarly, after each sub-step, = 

for level I + 1. Notice that the par¬ 
ticular sign depends on the entry in the mask, e.g. 
-l-x-surfaces refer to a positive sign in 0- 

3. When the parent and the child level are aligned in 
time, we sum up the contributions and correct the 
cell values with. 


£)(0 ^ £)(0 + 


Vi 


E 


Vi+i 


( 8 ) 


We observe that the effectiveness of the algorithm de¬ 
pends crucially on the specific RP operators. For hydro¬ 
dynamics fields the R step is conservative if the operation 
is performed using local averages, which are second order 
accurate, 0{h?). Similarly, a safe choice for the P step 
is linear interpolation using limiters in order to control 


TABLE IL Summary of the combinations for restriction 
(R), prolongation (P), and mass correction (C) used in this 
work. AVG indicates average, LAG Lagrangian interpolation, 
WENO, WENOZ the interpolation method of [80l[83]. The 
order of convergence is reported for each RP operation. 


Name 

R 

order 

P 

order G 

a2e2 

AVG 

2 

ENO 

2 

/ 

a2e2n 

AVG 

2 

ENO 

2 

/ 

a2wz6 

AVG 

2 

WENOZ 

6 

/ 

a2wz6n 

AVG 

2 

WENOZ 

6 

(Zf 

1414 

LAG 

4 

LAG 

4 

/ 

1414n 

LAG 

4 

LAG 

4 

/ 

w4w4 

WENO 

4 

WENO 

4 

/ 

w4w4n 

WENO 

4 

WENO 

4 

/ 


oscillations. However, for neutron star spacetime simula¬ 
tions high-order operators may be important for accuracy 
and faster convergence. As indicated in Tab. [ITj we have 
implemented several RP operators, including ENO 2nd 
order [82], Lagrangian and WENO 4th order [39]. In the 
next sections we will present results for various combina¬ 
tions of RPC operators. 


C. Atmosphere treatment 

For the simulation of neutron star spacetimes, the vac¬ 
uum region outside the stars requires special treatment. 
As described in [33], we use a low-density static and 
barotropic atmosphere at a density level 

Patm = fatm max[p(t = 0)]. (9) 

During the recovery of the primitive variables from the 
conservative variables, a point is set to atmosphere if the 
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density is below the threshold 


Pthr fthr Patm- 


( 10 ) 


The atmosphere treatment violates mass conservation 
and can potentially affect, or invalidate, the improve¬ 
ments related to the conservative AMR. In the follow¬ 
ing we will investigate this aspect in some detail ex¬ 
perimenting with parameters fatm G [10~^^,10 ~^] and 
fthr = (10^, 10^, 10^, 10'^) (see in particular Sec. V 31. 


IV. SIMULATION ANALYSIS 

In this section we describe the quantities employed for 
the analysis of our simulations. Notably, we introduce 
some diagnostic which are helpful to investigate the en¬ 
ergetics/geometry of the ejected material. 

The performances of the conservative AMR scheme are 
mainly tested using the baryonic, or rest-mass, mass in¬ 
tegral. 


Mb = J d^x = J d^x ^ D , (II) 

which should remain constant during the evolution, com¬ 
pare Eq. The rest-mass, and the other integrals dis¬ 
cussed in this work, are calculated on each refinement 
level. We usually report results for a given level, which is 
the appropriate one for the particular quantity; e.g. the 
baryonic mass is reported on the I = 1,2 level. 

The merger remnant of several BNS configurations 
considered here is a hypermassive neutron star (HMNS) 
which collapses to black hole on a dynamical timescale. 
The lifetime thmns is typically calculated from the mo¬ 
ment of merger (see below) to the time an apparent hori¬ 
zon forms. The black hole is then characterized by its 
horizon mass Mbh and spin jbh computed from the ap¬ 
parent horizon with average radius xah- 

The rest-mass of the accretion disk that forms after 
collapse is computed as. 


Mdisk = / d^x , (12) 

J i->rAH 

where the domain of integration excludes the spherical 
region inside the apparent horizon. 

The ejected material is defined by the two conditions, 

Ut < —1 and Vr = v^Xi > 0 , (13) 


been used in [JB] but not in [^. The total ejecta mass 
is computed as. 


Mejecta = / d^X , (14) 

Ju 

where the integral is computed on the region, 

U = {x* = {x,y,z) : Ut < —1 and Xj. > 0} , (15) 

on which material is unbound according to ITU). 

In order to investigate the energetics and geometry of 
the ejecta we consider different sets of integrals in the 
(x, ?/)-plane, in the (x, zj-plane, and the full 3D-domain. 
The kinetic energy of the ejecta can be approximated as 
the difference between the total energy Aejecta (excluding 
gravitational potential energy), and the rest-mass and 
the total internal energy C/ejecta [M] . 

^ejecta — E/ejecta (Afejecta fAjecta) 

= f d^x D{e — 1 — e) , (16) 

Ju 

where e = au^h — p/{pau*). Additionally, we compute 
the D-weighted integral of = ViV^, 

(17) 

(18) 

and the quantities. 


{v)p = 
{v)z = 


^ /^dxd^ 

\ fu D ) 




(P) 

{z) 


/ fi^dxdy D (x^ -by^) \ 

V D ) 

[ L dxdz n z^\ 

V ^ dxdz D J 


(19) 

( 20 ) 


(p) and (z) roughly estimate the geometric distribution of 
the ejecta. Similar integrals have been proposed in |54) . 
but in that case they were employed in three dimen¬ 
sions. We will use the approximation of the kinetic en¬ 
ergy Tejecta in Sec. VIIA and discuss the weighted veloci- 


ties {v)p^z and the (p), {z) for the case-study in Sec. VIII 
We compute the entropy “indicator' 


S = 


P 


( 21 ) 


where Ut = —W{a — j5iv‘^) is the first lower component of 
the fluid 4-velocity, and x* = {x,y,z). The first con¬ 
dition in (13) assumes fluid elements follow geodesics 
and requires that the orbit is unbound. This is a sim¬ 
ple criterion we use for continuity with previous work, 
e.g. [millj, and should at least capture the correct or¬ 
der of magnitude. The condition > 0 requires that the 
material has an outward pointing radial velocity; it has 


where T^ and Ki are locally determined by the value of 
p, see Tab. |llj In cases where the additional thermal con¬ 
tribution to the pressure Pth is small S' ^ 1, while in 
presence of shock heating S 1. 

Finally, gravitational waveforms are calculated via the 
curvature invariant 'I '4 and performing multipole decom¬ 
position on extraction spheres [32]. We work with the 
metric multipoles rhim, which are reconstructed from the 
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TABLE III: Grid and parameters configurations for single star 
tests. L denotes the total number of boxes, Z*"” is the finest 
non-moving level, n (n*"”) is the number of points in the 
fixed (moving) boxes, /iQ,hi_i are the grid spacing in level 
I — 0,L — 1 and , fatm the atmosphere level, and fthr the 
atmosphere threshold factor. The resolution in level I is hi = 
hol2K 


Single star test 

L 

jmv 

n 


ho 

hh-i 

fat 

m 

fthr 

TOVstatic 

5 

- 

56 

56 

2.0 

0.125 

10 ’ 

-9 

10 ^ 

TOVstatic 

5 

- 

56 

56 

2.0 

0.125 

10 “ 

11 

10 ^ 

TOVboost 

5 

- 

128 

128 

2.0 

0.125 

10 ’ 

-9 

10 ^ 

TOVboost 

5 

- 

128 

128 

2.0 

0.125 

10 “ 

11 

10 ^ 

TOV„.i9 

7 

- 

128 

128 

9.6 

0.150 

10 ’ 

10 

10 ^ 

TOV„.i9 

7 

- 

128 

128 

9.6 

0.150 

10 “ 

11 

10 ^ 

TOYmig 

7 

- 

128 

128 

9.6 

0.150 

10 ’ 

12 

10 ^ 

TOYmig 

7 

- 

128 

128 

9.6 

0.150 

10 ’ 

13 

10 ^ 

TOYmig 

7 

- 

128 

128 

9.6 

0.150 

10 “ 

10 

10 ^ 

TOYmig 

7 

- 

128 

128 

9.6 

0.150 

10 ’ 

11 

10 ^ 

TOVrnig 

7 

- 

128 

128 

9.6 

0.150 

10 - 

11 

10 ^ 

TOVmig 

7 

- 

128 

128 

9.6 

0.150 

10 “ 

11 

10 "^ 

RNSbut 

6 

1 

128 

64 

2.0 

0.0625 

10 ' 

-9 

10 ^ 

BUkcp 

7 

2 

144 

96 

4.0 

0.0625 

10 ’ 

-9 

10 ^ 


curvature multipoles using the frequency domain integra¬ 
tion of [84j . using the initial circular gravitational wave 
frequency as a cutting frequency, see Tab. [I^ All the 
waveforms are plotted against the retarded time, 

U = t-r^=t- Textr “ 2M In (rextr/2M - 1) , (22) 

where the extraction radius is fextr ~ 750 Mq. 


V. SINGLE NEUTRON STAR TESTS 




time t/M 0 



time t/M 0 

FIG. 3: Results of the TOVsta«c test. Top: Initial density 
profile of TOVstatic test along the a;-axis. The buffer zones of 
the refinement levels are shaded in gray. Middle: The relative 
rest-mass change 11— \ for different RPG combinations. 

Bottom: The time derivative of the rest-mass. 


Our conservative AMR implementation has been 
tested and validated in full-general relativistic simula¬ 
tions of single star spacetimes. In this section, we present 
five different tests, namely 


TOVstatic a static (nonrotating) neutron star with re¬ 
finement levels inside the star (Sec. VI); 

TOVfjoost a boosted, nonro tatin g neutron star crossing 
refinement levels (Sec. V 2); 


TOYmig a migration of an unstable spherical configu¬ 
ration to a stable one crossing refinement levels 
(Sec.[Vlt; 


For each test we perform simulations for different com¬ 
binations of the restriction (R), prolongation (P), and 
correction (C) step, as indicated in Tab. |nj The grid pa¬ 
rameters are reported in Tab. |III[ The most important 
quantity we are focusing on is the rest-mass. 

For the simulations we use both the BSSN and the Z4c 
system. Although no major differences between BSSN 
and Z4c are observed for what concerns mass conserva¬ 
tion, Z4c evolutions show overall smaller violations of 
Einstein constraints. 


1. TOVstatic 


RNSbut a uniformly rotating neutron star with refine¬ 
ment levels inside the star (Sec. V 4); 


RNSxep a neutron star close to the Kepler limit, which 
is perturbed and finally disrupted crossing refine¬ 
ment levels (Sec. V 5). 


We investigate a spherical star with a gravitational 
mass of 1.4M0. The initial data are calculated with a 
polytropic EOS with K — 100 and F = 2. The star 
is then evolved with the F-law EOS. The grid is pre¬ 
pared such that the finest refinement level I = A is fully 
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FIG. 4: Results of the TOVooost test. Top: Evolution of 
the density profile along the a:-axis; the profiles correspond 
to times t = 0, 50,100,150, 200, 250, SOOMq, and boost in the 
negative x direction. The buffer zones of the refinement levels 
are shaded in gray. Middle: The relative rest-mass change for 
different RPC combinations. Bottom: The time derivative of 
the rest-mass. 


contained in the star covering half diameter, and level 
I = 3 ends at the star surface. This is shown in the 
top panel of Fig. which collects the results. Although 
at the continuum the solution is trivial (static), in nu¬ 
merical simulations some dynamics is observed due to 
truncation errors. This is mostly triggered by the artifi¬ 
cial atmosphere treatment close to the star surface, and 
by truncation errors on the refinement levels I = 2,3. 
Thus, differences in the RPC steps influence the overall 
dynamics of the system. 

The middle and bottom panel of Fig. [^show the rel¬ 
ative error in the rest-mass and its time derivative. The 
conservative AMR (C step) improves the mass conserva¬ 
tion of about ~ 2 orders of magnitude, independently on 
the particular RP choice. Additionally, we observe differ¬ 
ences between the RP combinations. Even using C, the 
4th order WENO and Lagrangian RP introduce spurious 


FIG. 5: Results of the TOVmig test. Top: Evolution of the 
density profile along the a;-axis; the profiles correspond to 
times t = 0, 50,100,150, 200, 250, 3OOM0. The buffer zones 
of the refinement levels are shaded in gray. The star first 
expands reaching r ~ 50, then contracts, then bounces back 
and forth several times. The inset shows the time evolution 
of the central density. The vertical dashed lines refer to the 
times shown in this panel. Middle: The relative rest-mass 
change for different RPG combinations. Bottom: The time 
derivative of the rest-mass. 


oscillations in the rest-mass derivative (see green and or¬ 
ange solid lines). In general, using the average R leads 
to the smallest errors. 

These results refer to an atmosphere density Patm = 
10“®. We have experimented with the a2e2 RP setup 
and an atmosphere density of Patm = 10“^^. The result 
is shown as a black dashed line in the middle and bot¬ 
tom panel of Fig. A lower atmosphere significantly 
improves the mass conservation. In this test, the error 
in the rest-mass derivative related to the C step is about 
\dMi,/dt\ ^ 10“®, while the one related to the atmosphere 
treatment is about \dMb/dt\ ^ 10 “AtmPatm_ Hence, op¬ 
timal results can only be obtained with a proper combi¬ 
nation of RPC and {fatm, Patm)- 
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time t/MQ 

FIG. 6: Results of the TOVmig test: influence of the atmo¬ 
sphere parameters. In the legend, the first number represents 
fatm, the second number fthr- Solid lines correspond to simu¬ 
lations with a2e2 RPC; dotted lines to simulations with a2e2n 
RPC, i.e. without C step. 


2. TOVboost 

Initial data is prepared using the same star model as 
Sec. |V 1[ which is now boosted in the negative a;-axis 
direction^. The grid is prepared such that the star is 
initially entirely covered by the finest refinement level 
I = 4. During motion, the star crosses completely the two 
finest refinement levels, as shown in Fig.(top panel). 

As visible in Fig. (middle and bottom panels) the 
C step improves mass conservation in most of the cases, 
but here its effectiveness depends more significantly on 
the RP choice than in the TOVstatic test. In particular, 
the C step is not effective with WENO RP. The a2e2 and 
a2wz6 schemes perform best, indicating the importance 
of a conservative R. Similarly to the previous test, we test 
the role of the atmosphere parameters on the optimal 
a2e2 setup. Lowering the atmosphere by a factor 100 
improves mass conservation by a factor 10 in this case 
(see dotted black line). 


3. TOVmig 

We investigate an unstable single neutron star config¬ 
uration. Initially the central density is pc = 7.9934 • 10“^ 
and the gravitational mass 1.4476M0, see e.g. [39]. The 
star is in an unstable equilibrium, truncation errors trig¬ 
ger a migration to a stable configuration, which involves 
violent nonlinear oscillations on dynamical timescale. 
During these expansions and contractions, matter crosses 
the grid refinement levels. When matter reaches the grid 


^ We have further tested our implementation by boosting the star 
in all the directions, and both applying bitant symmetry, i.e. 
evolving only z > 0, and simulating the full numerical domain. 




time t/M0 



FIG. 7: Results of the RNSbu? test. Top: Density profile 
(red) and momentum density (blue) along the x-axis Middle: 
Relative rest-mass change for different RPC combinations. 
Bottom: The time derivative of the rest-mass. 


outer boundary some rest-mass falls out of the grid, but 
typically mass conservation is mostly affected by the in¬ 
teraction between the star low-densities outer layers and 
the atmosphere. Results are summarized in Fig. for 

Patm = 10-“ {fthr = 10^). 

We observe the conservative AMR is effective up to 
times t < AOQMq, that corresponds to ~ 2 bounces of 
the star core; up to the first bounce the rest-mass conser¬ 
vation improves of about two order of magnitude if the C 
step is used. At times t > 400 matter densities p ~ 10“® 
reach outer regions, where the resolution is dropped by 
a factor of 16 and interaction with atmosphere becomes 
significant. 

Figure [^summarizes our experiments with atmosphere 
parameters. Lowering pthr by an order of magnitude 
leads to an improvement of the mass-conservation by ap¬ 
proximately one order of magnitude for the beginning 
of the simulation, while for different Patm and the same 
Pthr the error stays the same, as expected. Relative 
rest-mass violation can be minimized up to 10“® using 
fatm = 10“^^ and fthr = 10^. One can notice that, if 
the C step is not applied and the atmosphere is small 
enough {patm ^ 10“^°), a dramatic mass violation hap- 
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FIG. 8: Results of the RNSkbp test. Left: Density evolution along directions x and 2 . Right: The relative rest-mass change 
and rest-mass time-derivative for different RPC combinations. 


pens as soon as matter crosses the first refinement bound¬ 
ary {t ~ IOOM 0 ), see dotted lines in the figure. The 
same does not happen with the C step. As time ad¬ 
vances, rest-mass conservation is progressively corrupted 
in all the cases due to the drop in resolutions in the outer 
region reached by the low-density star outer layers bounc¬ 
ing back and forth. 


4- RNSbu7 

Initial data is a stable uniformly rotating neutron star 
described by a polytropic EOS with K = 100 and T = 2, 
and with = 1.28 • 10“^, axes ratio 0.65, and gravi¬ 
tational mass 1.6655Mq, e.g. [5S]. The initial data are 
computed with the RNS code [Ml EH- The star is evolved 
with the T-law EOS, for about 6 periods. Results are 
shown in Fig. 

As in the previous tests, the C step improves the results 
in many cases; the best RP setup is a2e2. The 1414 and 
1414n RP perform equally good at late times. Surpris¬ 
ingly, the nonconservative w4w4n RP is here observed to 
give good results, and at the end of the simulation, it is 
comparable to a2e2. 


5. RNSkbp 

Initial data is a rotating neutron star at the Kepler 
limit modeled by a polytropic EOS with K = 100 and 
r = 2, and with pc = 1.444 • 10“^, axes ratio 0.58, and 
gravitational mass 1.7498M0. The star is evolved with 
the P-law EOS with P = 1.9; the lower polytropic ex¬ 
ponent triggers the star expansion with matter crossing 
several refinement levels. 

The left panels of Fig. [^ show how the matter expands 
along the x-axis and the z-axis over time, i.e. perpen¬ 
dicular and along the symmetry axis. The right panels 
of the figure show the mass conservation. The best RPC 
combinations are again a2wz6 and a2e2. 


6. Summary of single star tests 

Summarizing the results of the single star tests, we 
find the best mass conservation using the a2e2- scheme, 
i.e. the average restriction operation and a 2nd order 
ENO interpolation for the prolongation. The a2e2 sim¬ 
ulations show, on average, the smallest dMh/dt and no 
artificial oscillations in 1 — = 0). The lat¬ 

ter are present in at least one test for all other setups 
than a2e2. Additionally, the TOVstatic, TOVboost, and 
TOVmig tests suggest that also the artificial atmosphere 
treatment leads to mass violation. The stability of the 
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simulation improves with higher atmosphere values, but 
the mass conservation improves for lower atmosphere 
thresholds. An optimal setup is necessarily a compro¬ 
mise between these two effects. The largest violations 
of rest-mass conservations are observed in the lowest re¬ 
solved regions; where the violation becomes independent 
on the C step and the atmosphere values (i.e. it is mostly 
due to resolution). 


VI. BNS CONFIGURATIONS & GRID SETUP 


with resolution hi, and up to ri ^ rij2^~^ with resolu¬ 
tions hi = for I > 1. 

For all binary evolutions we run both the a2e2 and 
a2e2n RPC schemes. The a2e2 scheme is chosen because 
of its robustness and best performances in our previous 
tests; a2e2n is considered in order to assess the effect of 
the C step in the AMR strategy. We have not considered 
other combinations due to the computational overhead 
that they would imply. We set fatm = 10“^^ and fthr = 
10 ^ for all simulations. 


For this work we have prepared several BNS irrota- 
tional configurations in quasiequilibrium and circular or¬ 
bits; all the configurations are reported in Tab. IV Initial 
data are calculated with the LORENE [88] code. 

Our BNS sample spans the EOS sample of Tab. jn] 
and, for each EOS, two mass ratios^ q = Ma/Mb = 
1,1.16 are considered for a fixed total binary mass of M = 
Ma + Mb = 2.7. All EOS support maximal neutron star 
masses > 2Mq in agreement with recent observations [HU 
1 ^ . and the adiabatic speed of sound is Cg < c for a 
density range up to the maximum density supported by 
a stable TOV star. The compactnesses of the stars lie 
within C € [0.129,0.187]. The tidal coupling constant 
spans K 2 G [75,331] (see (23) for the definition). Notice 
that stiffer EOS have larger and, for the same EOS 
and M = 2.7, a larger q implies a larger The initial 
GW frequency of all binaries is MUJ 22 — 0.052. 

Additionally, we computed a, q = 1.5 and M = 2.5 con¬ 
figuration with the MSlb EOS (MSlb-100150). MSlb- 
100150 has a highly deformable EOS and = 461. The 
choice of parameters (EOS, q, M) of this configuration 
could be considered as “extreme” given the double pulsar 
population, e.g. m- However, the double pulsars sam¬ 
ple is rather small to be a significant statistics and the 
MSlb-100150 parameters are possible. 

Some of our M = 2.7 configurations have already been 
investigated in full general relativity in [S31[SU- Thus, 
the choice of initial data allows us to compare results 
with the literature. We will also compare with results 
of |93j employing smooth particles hydrodynamics and 
conformal flatness, although the evolution method differs 
and initial data are not prepared in the same way as here. 

For the BNS evolutions we use the Z4c scheme [67ll68| . 
and constraint preserving boundary conditions |68[ I94j . 
For all our runs a grid consisting of L = 7 refinement 
levels is used, levels with I > Z""' = 4 are dynami¬ 
cally moved. The grid spacing and outer boundary posi¬ 
tion depends on the employed model and is reported in 
Tab. jlVj For refinement level Z = 0 we employ the spheri¬ 
cal patches, as described in Sec. but we do not evolve 
matter on them. Indicating with r; the inradius on re¬ 
finement level I > 0, GRHD is evolved up to ri ~ n-/ii/2 


^ We define the mass-ratio to be always f? > 1. 


VII. BNS MERGERS WITH M = 2.7 AND 

g = 1,1.16 

We first discuss our results for configurations with mass 
ratios g = 1,1.16 and a total binary mass M = 2.7. We 
focus on the effect of mass ratio and the EOS on the 
merger dynamics, ejecta and gravitational waves. Also, 
we show the use of conservative AMR significantly im¬ 
proves the simulation of the merger remnant. Several 
results are reported in Tab. |V| and collected in Fig. |U 
Fig. [^ Fig. [^and Fig. [^ to which we refer during the 
discussion. 


A. Effect of EOS and q on merger dynamics 


The initial configurations, prepared in quasicircular or¬ 
bits at the same GW frequency MW 22 = 0.052, evolve 
for about 3 to 5 orbits before merger, depending on the 
EOS and mass ratio q. Here, the moment of merger is 
defined as the time corresponding to the peak of 
the £ = m = 2 multipole of the GW amplitude (see be¬ 
low). Stiffer EOSs give shorter inspiral (less revolutions) 
and lower dimensionless GW frequency at merger, see 
in Tab. 0 Unequal-mass systems are character¬ 
ized by slightly shorter inspiral than equal-mass ones and 
smaller merger frequencies of about ~ 3%. These prop¬ 
erties can be understood considering the values of the 
main (£ = 2) tidal polarizability parameter (tidal cou¬ 
pling constant hereafter) [M] , 


kI = 2 


f k2 
\{i+qrc% 


_^kl\ 
{l + qrc%j 


(23) 


where are the 1 = 2 dimensionless Love numbers 
of the individual stars [9611M] . in our sample. The re¬ 
sults agrees with the analysis of |100j . Essentially, for 
the same mass, stars with stiff EOS have larger radii 
than those with soft EOS, and attractive tidal interac¬ 
tions are stronger for larger values of ! thus, stiffer 
EOS binaries merge at lower frequencies. Notice that: 
(i) g > 1 configuration have slightly larger values of 
than q = 1] (ii) in our sample of configurations, EOS 
effects are typically larger than mass-ratio effects. The 
late-inspiral dynamics and GWs have been subject of re¬ 
cent work, e.g. mm and we do not discuss them any 
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TABLE IV: Initial BNS configurations and grid setup. First column defines the configuration name. Next 11 columns describe 
the physical properties: EOS, gravitational mass of the individual stars Ma,b, baryonic mass of the individual stars Mb a,b, 
stars’ compactnesses Ca,b, tidal coupling constant initial gravitational wave circular frequency MUJ 22 , ADM-Mass Madm, 
ADM-angular momentum Jadm- Next 8 columns describe the grid conhguration: finest grid spacing Hl-i, radial resolution 
inside the shells hr, number of points in the fix (moving) n (n™'”) levels, radial point number Ur and angular point number ng 
in the shells, inradius up to which GRHD equations are solved ri, and boundary rg. Notice that we divide the configurations 
in 3 different grid setups Rl, R2, R3 (compare simulation name). All conhgurations are evolved with and without the C step, 
which we denote with a “c” or “n” in the configuration name. 


Name 

EOS 

Ma 

Mb 

Mb A 

Mb B 

Ca Cb 

T 

f^2 

MuJ22 

Madm 

Jadm 

hh-i 

hr 

n 

n™” 

Ur 

ne 

ri 

rb 

MS1-135135-R2C 

MSI 

1.35 

1.35 

1.46 

1.46 

0.139 0.139 

325 

0.052 

2.676 

7.16 

0.240 

7.68 

160 

80 

160 

70 

614 

1870 

MSl-135135-R2n 

MSI 

1.35 

1.35 

1.46 

1.46 

0.139 0.139 

325 

0.052 

2.676 

7.16 

0.240 

7.68 

160 

80 

160 

70 

614 

1870 

MS1-125145-R2C 

MSI 

1.45 

1.35 

1.61 

1.38 

0.148 0.129 

331 

0.052 

2.673 

7.10 

0.240 

7.38 

160 

80 

160 

70 

590 

1870 

MSl-125145-R2n 

MSI 

1.45 

1.25 

1.61 

1.38 

0.148 0.129 

331 

0.052 

2.673 

7.10 

0.240 

7.38 

160 

80 

160 

70 

590 

1870 

H4-135135-R2C 

H4 

1.35 

1.35 

1.47 

1.47 

0.147 0.147 

210 

0.052 

2.674 

7.13 

0.2232 

7.1424 

160 

80 

160 

70 

571 

1739 

H4-135135-R2n 

H4 

1.35 

1.35 

1.47 

1.47 

0.147 0.147 

210 

0.052 

2.674 

7.13 

0.2232 

7.1424 

160 

80 

160 

70 

571 

1739 

H4-125145-R2C 

H4 

1.45 

1.25 

1.59 

1.35 

0.158 0.136 

212 

0.052 

2.674 

7.10 

0.230 

7.36 

160 

80 

160 

70 

589 

1792 

H4-125145-R2n 

H4 

1.45 

1.25 

1.59 

1.35 

0.158 0.136 

212 

0.052 

2.674 

7.10 

0.230 

7.36 

160 

80 

160 

70 

589 

1792 

ALF2-135135-R2C 

ALF2 

1.35 

1.35 

1.49 

1.49 

0.161 0.161 

138 

0.052 

2.675 

7.15 

0.202 

6.464 

160 

80 

160 

70 

517 

1574 

ALF2-135135-R2n 

ALF2 

1.35 

1.35 

1.49 

1.49 

0.161 0.161 

138 

0.052 

2.675 

7.15 

0.202 

6.464 

160 

80 

160 

70 

517 

1574 

ALF2-125145-R2C 

ALF2 

1.45 

1.25 

1.61 

1.37 

0.172 0.150 

140 

0.052 

2.673 

7.08 

0.200 

6.4 

160 

80 

160 

70 

512 

1558 

ALF2-125145-R2n 

ALF2 

1.45 

1.25 

1.64 

1.37 

0.172 0.150 

140 

0.052 

2.673 

7.08 

0.200 

6.4 

160 

80 

160 

70 

512 

1558 

SLy-135135-R2c 

Sly 

1.35 

1.35 

1.49 

1.49 

0.174 0.174 

74 

0.052 

2.675 

7.15 

0.1824 5.8368 

160 

80 

160 

70 

467 

1421 

SLy-135135-R2n 

Sly 

1.35 

1.35 

1.49 

1.49 

0.174 0.174 

74 

0.052 

2.675 

7.15 

0.1824 5.8368 

160 

80 

160 

70 

467 

1421 

SLy-125145-R2cl 

Sly 

1.45 

1.25 

1.62 

1.38 

0.187 0.161 

75 

0.052 

2.673 

7.07 

0.1824 5.8368 

160 

80 

160 

70 

467 

1421 

SLy-125145-R2nl 

Sly 

1.45 

1.25 

1.62 

1.37 

0.187 0.161 

75 

0.052 

2.673 

7.07 

0.1824 5.8368 

160 

80 

160 

70 

467 

1421 

SLy-125145-R2c2 

Sly 

1.45 

1.25 

1.62 

1.37 

0.187 0.161 

75 

0.052 

2.673 

7.07 

0.188 

6.106 

160 

80 

160 

70 

488 

1464 

SLy-125145-R2n2 

Sly 

1.45 

1.25 

1.62 

1.37 

0.187 0.161 

75 

0.052 

2.673 

7.07 

0.188 

6.106 

160 

80 

160 

70 

488 

1464 

MSlb-100150-Rlc 

MSlb 

1.50 

1.00 

1.64 

1.06 

0.157 0.109 

461 

0.042 

2.479 

6.16 

0.291 

9.312 

128 

64 

128 

56 

596 

1820 

MSlb-100150-Rln 

MSlb 

1.50 

1.00 

1.64 

1.06 

0.157 0.109 

461 

0.042 

2.479 

6.16 

0.291 

9.312 

128 

64 

128 

56 

596 

1820 

MSlb-100150-R2c 

MSlb 

1.50 

1.00 

1.64 

1.06 

0.157 0.109 

461 

0.042 

2.479 

6.16 

0.2328 

7.4496 

160 

80 

160 

70 

596 

1814 

MSlb-100150-R2n 

MSlb 

1.50 

1.00 

1.64 

1.06 

0.157 0.109 

461 

0.042 

2.479 

6.16 

0.2328 

7.4496 

160 

80 

160 

70 

596 

1814 

MSlb-100150-R3c 

MSlb 

1.50 

1.00 

1.64 

1.06 

0.157 0.109 

461 

0.042 

2.479 

6.16 

0.194 

6.208 

192 

96 

192 

84 

596 

1810 

MSlb-100150-R3n 

MSlb 

1.50 

1.00 

1.64 

1.06 

0.157 0.109 

461 

0.042 

2.479 

6.16 

0.194 

6.208 

192 

96 

192 

84 

596 

1810 


further here. In the following we focus on the postmerger 
dynamics. 

The postmerger dynamics has a rich phenomenology 
depending on the main binary properties: total mass, 
mass-ratio, EOS and stars’ spin (see e.g. [53l [Ml [56l 
isiiiM una for recent work). In the case of irrota- 
tional binaries and M = 2.70Mq, equal-mass mergers 
result in a massive differentially rotating compact ob¬ 
ject, which oscillates violently (see the Pmax = max(p) 
evolution in Fig. 10 right after merger). The compact 


object’s angular momentum is redistributed from the in¬ 
ner region to outer ones by torque and nonlinear hydro- 
dynamical interaction. The stability of the object de¬ 
pends on the mass, EOS and dissipative processes (see 
below). Following the literature |104j . we define this ob¬ 
ject as a hypermassive neutron star (HMNS), in case 
its rest-mass is larger than the maximum rest-mass of 
a stable uniformly rotating star described by the same 
EOS, or a supramassive neutron star (SMNS), in case 


its rest-mass is smaller. If the object does not exceed 
the rest mass of a stable TOV-solution, we simply refer 
to it as massive neutron star (MNS). These definitions 
apply to equilibrium configurations, in particular to cold 
EOS and axisymmetry, hence, although of common use, 
they cannot be rigorously applied to the merger rem¬ 
nants. In most cases HMNS are objects that are dynam¬ 
ically unstable and collapse to a black hole on timescales 
of ^ 2000 — IOOOOMq ~ 10 — 50 ms; whereas SMNSs are 
objects that appear stable on those timescales, but may 
eventually collapse later on due dissipative processes, 
e.g. loss of angular momentum radiated via GWs. Snap¬ 
shots of the density distribution and velocities in the or¬ 
bital plane are presented in Fig.[^ the simulation time is 
close to the moment of merger. 

Three of our <7=1 configurations, H4-135135, ALF2- 
135135, and SLy-135135, merge in a HMNS which col¬ 
lapses to a black hole (BH) within thmns 2000—5000 ^ 
10 — 25 ms from the merger moment. All these EOSs sup- 
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TABLE V: Summary of the numerical results for the M — 2.7M0-simulations. Columns: Simulation name, merger time, 
merger frequency (stated dimensionless and in Hz), final remnant, the lifetime of the HMNS thmns stated in solar masses and 
milliseconds, 2nd peak fs- and /2-mode frequency (dimensionless and in Hz), mass and kinetic energy of the ejected material 
Mejecta (see Fig. 10 1 , the mass of the disk surrounding the central object Mdisk measured ~ 2OOM0 after BH formation, the 
black hole mass ALbh and its dimensionless angular momentum /bh- 


Name 

tmrg 


Remnant thmns Mlj 22 


MuJ22 

/2 

-^^ejecta 

^ejecta 

Mdisk Mbh 

JBH 




[kHz] 

(ms) 

[kHz] 


[kHz] 

[10-3] 

[10"^] (10®° erg) 

[10-°] 



MSl-135135-R2c 

1479 

0.112 

1.38 

MNS - 0.134 

1.60 

0.166 

1.99 

0.7 

0.1 (0.2) 

- 

- 


MSl-135135-R2n 

1476 

0.114 

1.36 

MNS - 0.135 

1.61 

0.170 

2.04 

1.2 

0.1 (0.2) 

- 

- 


MS1-125145-R2C 

1420 

0.110 

1.32 

MNS - 0.130 

1.56 

0.172 

2.06 

5.8 

0.7 (1.2) 

- 

- 


MSl-125145-R2n 

1419 

0.111 

1.33 

MNS - 0.125 

1.50 

0.157 

1.88 

3.2 

0.2 (0.4) 

- 

- 


H4-135135-R2C 

1804 

0.129 

1.54 

HMNS-^BH 5130 (25) 0.146 

1.75 

0.214 

2.57 

0.6 

0.3 (0.5) 

10.8 

2.48 

0.62 

H4-135135-R2n 

1803 

0.130 

1.55 

HMNS-^BH 4470 (22) 0.145 

1.73 

0.216 

2.58 

0.6 

0.3 (0.6) 

8.5 

2.54 0.65 

H4-125145-R2C 

1822 

0.120 

1.44 

HMNS - 0.140 

1.68 

0.197 

2.36 

6.0 

1.6 (2.8) 

- 

- 

- 

H4-125145-R2n 

1820 

0.120 

1.44 

HMNS - 0.146 

1.75 

0.194 

2.32 

4.0 

1.2 (2.3) 

- 

- 

- 

ALF2-135135-R2C 

2148 

0.142 

1.71 

HMNS^BH 3760 (19) 0.168 

2.01 

0.235 

2.81 

3.5 

0.4 (0.7) 

17.8 

2.43 

0.62 

ALF2-135135-R2n 

2145 

0.142 

1.71 

HMNS-^BH 3770 (19) 0.165 

1.98 

0.230 

2.75 

2.0 

0.4 (0.7) 

21.1 

2.44 0.63 

ALF2-125145-R2C 

2028 

0.138 

1.65 

HMNS - 0.157 

1.88 

0.222 

2.66 

3.9 

0.4 (0.8) 

- 

- 

- 

ALF2-125145-R2n 

2027 

0.139 

1.66 

HMNS - 0.160 

1.91 

0.225 

2.69 

10.6 

1.0 (1.9) 

- 

- 

- 

SLy-135135-R2c 

2504 

0.168 

2.01 

HMNS-^BH 2159 (11) 0.206 

2.46 

0.292 

3.49 

12.2 

4.0 (7.1) 

8.4 

2.48 

0.64 

SLy-135135-R2n 

2495 

0.168 

2.01 

HMNS^BH 2577 (13) 0.207 

2.48 

0.290 

3.47 

14.2 

5.9 (10.5) 

9.6 

2.49 

0.64 

SLy-125145-R2cl 

2353 

0.162 

1.93 

HMNS-^BH 3020 (15) 0.184 

2.20 

0.286 

3.42 

6.5 

2.8 (5.1) 

17.9 

2.40 

0.58 

SLy-125145-R2nl 

2350 

0.161 

1.93 

HMNS^BH 2870 (14) 0.187 

2.24 

0.283 

3.39 

4.5 

1.7 (3.0) 

14.5 

2.46 

0.61 

SLy-125145-R2c2 

2350 

0.161 

1.92 

HMNS^BH3310 (16) 0.186 

2.23 

0.285 

3.41 

6.2 

2.1 (3.7) 

18.4 

2.40 

0.58 

SLy-125145-R2n2 

2348 

0.160 

1.91 

HMNS-^BH2180 (11) 0.184 

2.20 

0.283 

3.39 

5.4 

2.5 (4.5) 

11.1 

2.49 

0.62 


port approximately the same maximum mass regarding 
single spherical stars, but the stiffer the EOS, the longer 
is Thmns- This fact can be understood by the following 
considerations. The range for the tidal coupling constant 
is ^2 £ [75,331], where soft (stiff) EOS binaries corre¬ 
spond to small (large) values in this range. Stiff EOS 
binaries are gravitationally less bound systems than soft 
EOS binaries: their binding energy at merger is larger 
(less negative) as well as the angular momentum. As a 
result, the HMNS has more angular momentum support 
at formation. However, the initial angular momentum is 
not the only factor that determines the lifetime of the 
HMNS. At formation, the HMNS density in the star core 
increases, the pressure response depends on the (effec¬ 
tive) adiabatic index of the fluid which is different for 
each EOS. As a result, the HMNS nonlinear oscillations 
and the efhciency of the angular momentum redistribu¬ 
tion depend on the EOS [52] • Stiffer EOSs have larger 
pressure support against gravity, especially at high den¬ 
sities. Finally, in a more realistic situation than the one 
simulated here (and on longer timescales), thermal sup¬ 
port, angular momentum transport driven by magnetic 
fields^ and cooling mechanisms (neutrinos) are expected 


® We notice the largest simulations with present techniques and 


to play a role. The lifetimes of the HMNS are stated in 
Tab. [Vj and our results agree with [52| within ±5 ms. 

The merger of MSl-135135, differently from the other 
q = 1 configurations, produces a differentially rotating 
object that is stable over the whole simulation time, 
i.e. 6000Mq ~ 30 ms after merger. Non-rotating stars 
described by the MSI EOS can support a maximum rest- 
mass of ^ 2.767M0. According to the previous defini¬ 
tion, we classify the merger remnant for the MSI models 
as a MNS. Considering the physics simulated here, we 
expect that the merger remnant will stabilize via GW 
emission reaching a uniformly rotating and cold config¬ 
uration on the characteristic timescale, tqw ~ J/T ~ 
{R)^/{Mf « 4 OOOOM 0 « 200 ms. 

The unequal-mass q = 1.16 configurations H4-125145 
and ALF2-125145 have a different merger remnant than 
the corresponding q = 1 configurations. In these cases we 
find an object stable over 5 OOOM 0 ^ 25 ms, but since the 
mass is still larger than the supported mass of the uni¬ 
form rotating model, remnants are HMNSs. We expect 
these configurations will collapse within t < tgw- Ref¬ 
erences IMIE2 found that similar configurations with a 
slightly different thermal component Tth = 1.8 form BHs. 


resolutions have not properly resolved magnetic field amplifica¬ 
tion effects [105]. 
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1-3 


-5 





SLy-135t 34R2C - x-y-plane t=2524.4 



-5 


-6 




SLy-125145R2c1 - x-y-plane 1=2393.1 


FIG. 9: 2D snapshot of density and velocity on the orbital plane shortly after the moment of merger. The velocity pattern 
is indicated by red arrows. The region inside the black contours contain unbound material on a logarithmic scale with 
pejecta = (10“, 10“®, 10“®, 10“^, 10“®, 10“®). The logarithm of the density log^Q (p) is visualized according to the color 
bar. Left (from top to bottom): MSl-135135-R2c, H4-135135-R2c, ALF2-135135-R2c, SLy-135135-R2c. Right (from top to 
bottom): MSl-125145-R2c, H4-125145-R2c, ALF2-125145-R2c, SLy-125145-R2cl. 
















































































































15 



FIG. 10: Evolution of several dynamical quantities for M = 2.7Mq q = 1, 1.16 configurations. Results for different EOS 
are in different color. For each configuration, the panel contains four plots. From top to bottom: rest-mass violation SM = 
Mb{t) — Mfjit = 0) on level f = 1; maximum density Pmax = max(p) on the grid scaled to its initial value Pmax(i)/pmax(t = 0); 
rest-mass of the ejected material Mejscta', kinetic energy of the ejecta Tsjecta- Results for the conservative AMR are presented 
with solid lines, while the corresponding results for the nonconservative AMR are shown with dashed lines. Vertical lines 
represent the moment of merger, i.e. tmrg determined by the maximum in \rh 22 \- 


A similar dynamics as for the MSl-135135 is observed in 
the merger of MSl-125145, where a stable MNS is pro¬ 
duced. The SLy-125145 forms, as in the q = 1 case, a 
black hole, but, following the general trend, the HMNS 
lifetime is longer. 

Due to the unequal mass ratio, the merger remnant is 
typically more deformed than the corresponding q = I 
and strongly non-axisymmetric at formation, see Fig. 
Unequal-mass binaries have more stable merger remnants 
than corresponding equal-mass ones (e.g. larger thmns)- 
The q — 1.16 HMNS/MNS are characterized by slightly 
larger radii than the q = 1 ones, and a different cen¬ 
tral density, Fig. Additionally, the mass-ratio has an 
effect on the ejecta as we shall see below. 


At formation, all the merger remnants show violent 
oscillations, visible in the evolution of Pmax in Fig. 10 


The softer the EOS, the larger are the oscillations, see in 
particular the SLy panels in the hgure. This property re¬ 
flects the pressure response of the EOS for density jumps 


around p> P 2 (Of. above and also [52]). The oscillations 
have a quasi-radial character, and relax either within few 
radial periods or before the onset of collapse. 

In cases with BH formation, the BH masses are of or¬ 
der 2.4 —2.5 Mq, and the dimensionless BH spin is of the 
order 0.58 — 0.64 for all the configurations. The evolu¬ 
tion of the BH parameters is presented in Eig. 11 (top 
and middle panels). These results suggest that, in this 
scenario, the BH formation and properties are mostly 
determined by the total mass of the system and depend 
only weakly on other details. However, uncertainties on 
these numbers are of the order of ~ 2% — 5%, and it is 
difficult to draw precise conclusions. 

The final BH is surrounded by an accretion disk of 
rest-mass M^isk ~ 0.05 — 0.2 Mq, see Tab. |V| and the 
bottom panel of Eigure 11 The disk geometry is es¬ 
sentially axisymmetric for all the configurations. During 
the evolution, the maximum density inside the disk de¬ 
creases from ^ 10“^ to ~ 10“^. At the moments the 
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FIG. 11: Black hole and disk evolution for simulations with 
and without conservative AMR. Top: black hole horizon mass. 
Middle: black hole dimensionless angular momentum. Bot¬ 
tom: disk rest-mass. 


BH masses and spins reach their plateaus (late times in 
our simulations), the dense regions of the disk extend 
up to distances < 30 ~ 45 km. Lower density, gravi¬ 
tationally bound regions larger than patm extend up to 
^ 100 — 130 ~ 150 — 200km. The accretion rate is of the 
order Mdisk 10 “®. 


B. Assessment of conservative AMR 

Before continuing the analysis of the physical prop¬ 
erties of the merger remnant we discuss here the accu¬ 
racy improvements due to the numerical algorithm de¬ 
scribed in Sec. |IIIB| Figure reports result obtained 
with (solid lines) and without (dashed lines) the C step 
in the AMR algorithm. The information of the figure 
is complemented with the entries of Tab. |V] For all the 
configurations the C step is crucial for the simulation ac¬ 
curacy after merger. 

Let us first discuss rest-mass conservation. As pointed 
out in the introduction, the rest-mass can in general in¬ 
crease or decrease. In our BNS simulations we identify 
two main and competitive causes for the violation of con¬ 
servation: (i) when fluid crosses refinement boundaries 
rest-mass tends to increase, (ii) the artificial atmosphere 
treatment tends to decrease the rest-mass. Clearly, the 
C step can improve only violations of type (i). 

For most of the configurations the use of the C step 
leads to an improvement of a factor of ~ 5, except for 
the MS1-135135-R2 configuration where an improvement 


by more than a factor of ^ 10 is observed. The only 
simulation were no significant improvement is observed 
is SLy-125145-R2, where the violation is < 20% from 
merger to the end of the run. 

Overall, the data show some dependence on the EOS. 
Without C step the mass conservation is in general bet¬ 
ter for softer EOS; this is probably related to the smaller 
star deformations. On the contrary, with C step, slightly 
larger errors are observed for softer EOS. We suggest 
that this is caused by the influence of numerical viscos¬ 
ity, that, in these runs, is more significant than in the 
runs without C step due to better overall conservation. 
Notice that the performance of the conservative AMR 
algorithm is always better than (or at most comparable 
to) the corresponding simulations without C step. In [53] 
we have employed larger grid boxes without C step in an 
attempt to optimize the performances of the nonconser¬ 
vative AMR for the remnant simulation. In Appendix [^ 
we present some experiments along this line showing that 
conservative AMR is, in general, a better strategy. 


Mass-violations influence the behavior and lifetime of 
the merger remnant, as evident from Fig. |10[ We observe 
systematic shifts in the collapse time of several HMNS 
although there are no qualitative differences due to the 
sufficiently high grid resolutions of our runs. For H4- 
125145 the mass violation in the outer layers in the H4- 
125145-R2n run determines a slightly different evolution 
of the MNS and a lower Pmax- 

We observe maximal differences of a factor of 3 in the 
ejecta mass measured on the coarsest level {I = 1 ) be¬ 
tween the runs with and without the C step. Figure [TOj 
(bottom panel for each EOS) shows that the differences is 
larger either shortly after merger time or at later times: 
no clear trend is identifiable. Thus, low density ejecta 
remain challenging to simulate even with conservative 
AMR (as long as nested boxes are used as opposed to lo¬ 
cal AMR tracking the ejecta). In particular, the artificial 
atmosphere influence is probably significant: (i) during 
inspiral we observe some spurious ejecta due to atmo¬ 
sphere fluctuation, and (ii) at late times, when ejecta 
have expanded into larger radii (coarser resolutions) we 
expect an effect as the one discussed for the TOV^ig test 
in Sec. IV 31 

Differences in the black hole and disk remnant are also 
observed, see Tab. jV] and Fig. 11 If the C step is not 
employed the estimated disk mass M^isk changes up to 
~ O.O 6 M 0 . In all configurations the final black hole mass 
and spin is overestimated when no C step is applied, 
which is probably related to the increase of the rest-mass 


visible in the upper panels (for each EOS panel) of Fig. 10 
(dashed lines). 

Finally, we mention that the GWs calculation during 
the inspiral is basically not influenced by the use of the C 
step. This is due to the fact that we have not attempted 
to refine the grid inside the star during that phase. Dur¬ 
ing orbital motion the stars stay compact and there is 
no need of further improving mass conservation. GWs 
in the post merger reflect the slightly different dynamics. 
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FIG. 12: Gravitational waves signals for M = 2.7 g = 1,1.16 configurations. For each configuration, the panel contains two 
plots. Top: K(rh,22); Bottom: Mu 22 - Results for the conservative AMR are presented in solid lines, while the corresponding 
results for the nonconservative AMR are in dashed lines. Vertical lines mark the moment of merger, i.e. tmrg determined by 
the maximum in |r/i22|- 


but the characteristic frequencies (see below) are essen¬ 
tially unaffected. 

In the following we will discuss exclusively the results 
obtained with conservative AMR scheme. 


C. Ejecta 

In this section we discuss the EOS and mass-ratio ef¬ 
fect on the dynamical ejecta. A detailed analysis of the 
dynamical formation of the ejecta will be presented in 
Sec, 

Figure shows the evolution of the ejecta mass for 
the various configurations; Tab. [V| reports the maximum 
value. Ejecta peaks happen during and shortly after 
the merger moment; the ejecta rest-masses at this time 
are about Mejecta ~ 10 “^ Mq, and in some cases reach 

Mejecta ^ 10 ^ AIq. 

The amount of ejected material depends on the EOS 
and on the mass ratio. If g = I larger ejecta are observed 
for softer EOS. For a given EOS (but except for SLy 


EOS), g = I.I6 configurations have larger ejecta than g = 
I ones. Similarly, the kinetic energy estimate computed 
according to Eq. (16) is larger for softer EOS and larger 
g than for stiff ones. Our results for MSI, H4, and ALF2 
configurations essentially agree with 


We stress that ejecta computations are challenging. At 
present, mass conservation and artificial atmosphere are 
the main factors limiting the accuracy. This is evident in 
the case of SLy configurations. The results in Sec. |VII B| 
suggest that the evolution of this soft EOS is less accurate 
than the others (see also discussion in [M]). We believe 
this is the reason why Mejecta is larger for SLy-I35I35 
than for SLY-I25I45. The poor mass conservation in 
SLY-I25I45 certainly affects the ejecta calculation. No¬ 
tice also that a similar setup as SLy-I35I35 has been 
evolved in [S5]; there, the ejecta mass was estimated to 
lie in the range between (4 • 10“^, 6.4 • 10“^). 
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— MS1-135135 — ALF2-135135 

— H4-135135 — SLy-135135 




— MS1b-100150-R3 — ALF2-125145 

— MS1-125145 — SLy-125145 

— H4-125145 



0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

frequency [kHz] 


FIG. 13: Gravitational waves spectra for M = 2.7 q = 1 
(top) and q = 1.16 (bottom) configurations. The colors cor¬ 
respond to Fig. 12 MSlb (red), ALF2 (orange), H4 (green), 
SLy (blue). The thick lines refer to the entire GW-signal, 
while the thin lines include only the GWs emitted after the 
moment merger t > tmrg. Important frequencies are marked 
in the plot: /mrg (triangles), fs (circles), and /2 (diamonds). 
Additionally the MSlb-100150-R3c is added in the bottom 
panel. For this setup no clear /^-frequency is visible. 


D. Gravitational waves 


The £ = m = 2 multipoles of the GWs are shown in 
Fig. 12 for all the configurations. For each EOS, each 


panel shows the real part of the wave (top) and the in¬ 
stantaneous GW frequency (bottom). The vertical line 
in each panel marks the moment of merger, defined as 
the peak of the amplitude |r/i 22 |. 

The emission from the orbital motion is the char¬ 
acteristic chirping signal, in which frequency and am¬ 
plitude monotonically increase. At these separations, 
the dynamics is strongly affected by tidal interactions 
(parametrized by «:|’), and the GWs phase carry infor¬ 
mation about the EOS. A detailed and accurate semi- 
analytical modeling of the inspiral up to merger has been 
given in m- The chirp signal ends at the amplitude 
peak. 

After the merger moment, the amplitude instanta¬ 
neously drops down since the two stars merge in a sin¬ 


gle body which has, for one instant, a quasispherical ge¬ 
ometry [39] (see also the frequency spikes). The post¬ 
merger signal is mainly characterized by the nonlinear 
oscillations of the merger remnant. As discussed above 
and elsewhere, e.g. [M], the merger remnant can be ap¬ 
proximated by a compact star oscillating nonlinearly at 
the proper frequencies. The m = 2 /-mode with fre¬ 
quency /2 is the most efficient emitter of GWs, and it 
is strongly excited at formation. Thus, the GW emit¬ 
ted by the HMNS/MNS is dominated by this frequency. 
Looking at the frequency in Fig.[^ large oscillations are 
present right after the merger moment and correspond to 
the very nonlinear phase described in Sec. jVII A[ softer 
EOS show larger oscillations. During early stages of 
the HMNS/MNS evolution, different modes are excited, 
see e.g. the spectrogram in |53j . Nonlinearity results in 
mode couplings, the main ones being the combination 
f± = F ± f 2 between the quasiradial mode F and the 
/2 |106j . In cases where a MNS is formed (MSI EOS), 
the frequency oscillations relax quickly; the power in the 
f± channels decreases, and the frequency essentially set¬ 
tles on the /2 mode. In cases where a HMNS is formed, 
the frequency monotonically increases as a result of the 
star contraction prior to collapse. 

Let us finally discuss the GW spectra shown in Fig.fl^ 
The figure includes, for each configuration, the spectrum 
of the entire signal as a thick line and the spectrum con¬ 
sidering only the signal for t > f„irg as a thin line. Some 
of the relevant frequencies are marked with bullets: the 
frequency at the waveform amplitude peak /„ii.g (trian¬ 
gles), a frequency fs related to a secondary postmerger 
peak (circles), and the /2 frequency corresponding to the 
main postmerger peak (diamonds). Recently, there has 
been intense research about the identihcation and char¬ 
acterization of this postmerger GW spectrum frequen¬ 
cies [Ml |92l I106II110| . Eor most of the configurations, 
the /2 frequency is clearly identifiable. Note however the 
double peak for the MSI models. 

The /2 frequency is smaller for stiffer EOSs; for hxed 
EOS, q = 1.16 configurations have slightly smaller /2 
than q = 1. Our /2 values agree with [Ml[Ml[TTT] . 

The origin of the secondary peak is not well under¬ 
stood. fs appears mostly related to the very late inspi¬ 
ral phase: several fs peaks are not present, or strongly 
suppressed, if the PSD is computed using only times 
t > tmrg- However, for configuration SLyl35135 and 
H4135135 one can notice a clear secondary peak also in 
the PSD of the signal at times t > tmrg- We observe that 
the fs peaks generated by signals at times t > are 
suppressed for unequal-mass configurations {q > 1). Our 
values of fs are in good agreement with the frequencies 
called /i in [lllj . Our PSD analysis might be compat¬ 
ible with the interpretation of |110j according to which 
the peak of the spectrum close to fs is due to two dif¬ 
ferent effects: the nonlinear mode coupling /_ (that can 
be extracted clearly using the t > tmrg signal only), and 
motion of spiral arms during the last stage of the merger 
process (but mostly at times t < tmrg) at a frequency 
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TABLE VI: Summary of the numerical results for the MSlb- 
100150 simulation. Columns: Grid identifier, time at merger 
tmrg, GW frequency at merger stated dimensionless and in 
kHz, the peak frequency of the GW spectrum during the 
HMNS phase /2 stated dimensionless and in kHz, and maxi¬ 
mum mass of the ejected material Mejecta- 


Resolution 

^mrg 

[Mq] 


/mrg 

[kHz] 

MuJ22 

h 

[kHz] 

-^^ejecta 

[lO-^Mol 

Rlc 

2675 

0.086 

1.11 

0.137 

1.77 

32.6 

Rln 

2640 

0.085 

1.10 

0.139 

1.79 

27.8 

R2c 

2710 

0.086 

1.11 

0.141 

1.82 

27.7 

R2n 

2701 

0.085 

1.10 

0.140 

1.81 

29.4 

R3c 

2754 

0.088 

1.13 

0.145 

1.87 

29.9 

R3n 

2757 

0.088 

1.14 

0.142 

1.83 

28.3 


called there fspiral • 


VIII. THE MSlB-100150 CONFIGURATION 

In this section we consider the evolution of a configura¬ 
tion described by the MSlb EOS |112LI113] with q = 1.5 
and binary mass M = 2.5Mq. The individual stars have 
masses = 1.5 Mq andM^ = 1.0 Mq. This configura¬ 
tion has (to our knowledge) the largest mass-ratio simu¬ 
lated in numerical relativity. A g = 1.5 has been already 
simulated in [56] for the soft EOS APR, but no gravi¬ 
tational wave signal was computed. The specific MSlb- 
100150 configuration considered here has never been sim¬ 
ulated before. We focus on this case study to discuss in 
some detail the dynamical mechanism that generates the 
ejecta in the strong field region and the ejecta geometry 
at their formation. Eurthermore, we point out that in the 
BNS parameter space the combination of stiff EOS and 
large mass-ratio (and a system with a low mass ^ 1 Mq 
companion) produces a rather peculiar merger remnant 
in which a MNS is surrounded by a massive and extended 
accretion disk. 


A. Dynamics and Merger remnant 

Eigure shows a 3D rendering of the density 
p during the merger process at selected times t ^ 
2560,2957, 3200, 5440. Both the bound and unbound 
parts are shown, using an inverse color scale: from yel¬ 
low to light blue (bound p) and from blue to red (un¬ 
bound p). About 1.5 orbit before the moment of merger 
the stars come in contact; the companion {Mb = ^Mq) 
is very deformed by the tidal field of the primary star 
{Ma = 1.5 Mq). We observe the first mass ejection 
from the low-density outer layers of the companion, 
p ^ ~ 10® g cm“^ around this time (see the 

green/blue tail in the top left panel). At later times, 
the companion is partially disrupted: some material is 


captured into the primary and forms a hot and differen¬ 
tially rotating core; other material forms a tidal tail, see 
the top-right and bottom-left panels. Low density mate¬ 
rial p < 10“^ in the outer part of the tidal tail becomes 
unbound, and it is ejected from these regions during two 
main episodes. The higher density material, closer to the 
primary star, expands by centrifugal forces but remains 
bound. The final merger remnant is composed of a high 
density hot core surrounded by a thick accretion disk 
of rest-mass ~ 0.3 Mq and of radius ~ 35Mq « 55km 
(bottom-right panel). The remnant is not expected to 
collapse since the total binary rest-mass is smaller than 
the maximum rest-mass supported by this EOS for spher¬ 
ical configurations (Tab. E- 

The rest-mass of the total ejected material is about 
Alejecta ~ 0.03 M©. The large amount of mass ejected 
by this configuration offers the possibility to study with 
enough accuracy the ejecta formation process. 

We have checked our results against resolution consid¬ 
ering three different grid setups (Tab. VI) and excluding 
the C step in the AMR algorithm. In Eig. [^we present 
the mass conservation and the maximum density evolu¬ 
tions for all setups. The conservative AMR improves re¬ 
sults: by the end of the simulation and in the worse case, 
Mb is conserved up to 0.3% (1.7%) if the C step is (is 
not) applied. Larger differences in the Mb are observed 
among different resolutions for the nonconservative AMR 
runs than for the conservative AMR ones. Interestingly, 
the central density of the remnant is denser without the 
C step (bottom panel, compare previous section). Ab¬ 
solute uncertainties in the rest-mass conservation are of 
order 2.5 • 10“^ Mq by the end of the simulation, and are 
about a factor 10 smaller of Mejecta- 


B. Ejecta formation 


Let us discuss the dynamical process at the origin of 
mass ejection. We identify two main hydrodynamical 
mechanisms: (i) the torque exerted by the central two- 
cores structure on the tidal tail; and (ii) shock waves 
generated in the region between the two cores. Most of 
the unbound mass is ejected at times close to the moment 
of merger 2650 and around the orbital plane with 

a small opening angle of < 15°. From the first three pan¬ 
els of Fig. one can clearly observe that mass is ejected 
mostly from the tidal tails primarily of the companion 
star; the torque mechanism (i) is the dominant one. 

In order to further investigate mass ejection, we con¬ 
sider 2D plots of the rest-mass density p, velocity u*, 
and entropy indicator S, on the orbital {x, y,z = 0) 
and perpendicular {x,y = 0,z) planes. Fig. 16 The 
color map refers to logj^g S, white contour lines refer to 
p = (10“^, 10“®, 10“^, 10“"^, 10“^), arrows to the veloc¬ 
ity pattern, and regions delimited by black solid lines 
highlight unbound material with contour densities p = 
(10“^®, 10“®, 10“®, 10“^, 10“®) on a logarithmic scale. At 
time t ^ 1900, the revolution/rotation of the cores exerts 
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FIG. 14: The strong-field merger dynamics of MSlb-100150. The figure shows four snapshots of the bound and unbound 
density p at t ~ 2560 (top-left), t ~ 2957 (top-right), t ~ 3200 (bottom-left), and t ~ 5440 (bottom-right). All subplots contain 
the same contour range and the same part of the computational domain. The bound density p is shown on a logarithmic scale 
from 10“® (yellow) to 10“® (blue), and highlighted with contours for p = (10“®, 10“^, 10“®). The unbound material is shown 
on a logarithmic scale from 10“® (blue) to 10“® (red). Top-left: About 1.5 orbit before the moment of merger the stars come 
in contact. At t ~ 2560 the companion {Mb = 1 Mq, left) is deformed by the tidal field of the primary {Ma ~ 1.5M©, right). 
Ejecta originate from the tidal tail of the companion, and are emitted around the orbital plane. Top-right: At t ~ 2957, shortly 
after the moment of merger, the companion is already partially disrupted, most of the ejecta is emitted around this time. 
Bottom-left: At t ~ 3200 material is also ejected by the shock-heating-driven mechanism described in the text in a direction 
perpendicular to the orbital plane. On larger scales (not shown in the plot) ejecta appear anisotropically distributed around 
the orbital plane with an opening angle ~ 10°. Bottom, right: The merger remnant is composed of a MNS with a high density 
core surrounded by an accretion disk of rest-mass ~ O.SMq. The entire disk has a radius of ~ 35Mq « 55km. 


torque on the low-density outer layers of the companion 
star. This material gains enough energy to become un¬ 
bound and the ejection process stars. The ejected mate¬ 
rial expands with initial velocities {v)p ^ 0.3 and decom¬ 
presses. At this times also minor ejecta due to shocks 
occur (Fig. top left). Between t ^ tmrg ~ 2650 and 
t ~ 2900 mass is also ejected from the tidal tail of the 
primary star. The entropy has a spiral-like pattern in S 
(Fig. 161; the influence of the thermal pressure compo¬ 
nent Pth is larger in less dense regions. At t ~ 3000 we 
observe another significant event that causes mass ejec¬ 
tion. As clear from the bottom panels of Fig.[^ in this 
case the ejection is triggered by the shock wave generated 
between the two density maxima of the MNS. The fluid 
is heated up and driven outward by the thermal pressure 
(corresponding high entropy regions). The mass is ini¬ 
tially ejected in a direction roughly perpendicular to the 
orbital plane, but it falls back on the orbital plane and 
acquires angular momentum by torque. 


Figure [T^ quantifies mass, kinetic energy, and geom¬ 
etry of the ejecta. The rest-mass of the total ejected 
material is about Mejecta ~ 0.03 M©. Notice that, consis¬ 


tently with the discussion in previous sections, the mass 
decrease is mostly a numerical effect due both to reso¬ 
lution and atmosphere setup. The kinetic energy of the 
ejecta is Tejecta 3.2 • 10““^ ^ 2.9 • 10®°erg. Regarding the 
geometry, lower panel of Fig. we observe that mass ex¬ 
pands inside the orbital plane more rapidly than perpen¬ 
dicular to it. The analysis of the {p) and (z) curves sug¬ 
gests that the ejecta extends mainly around the equato¬ 
rial plane with an opening angle of d ~ arctanz/p ~ 10° 
(compare Eq.(19l and (20)). On large spatial scales, the 
geometry is anisotropic. 

The basic mechanisms (i) and (ii) identified in this case 
study are rather general and at the origin of mass ejec¬ 
tion also in other configurations. Thus, the geometrical 
and kinematic properties of dynamical ejecta at their for¬ 
mation described here are expected to be representative, 
at least at a qualitative level (see also [Ml El]). 

Clearly, configuration details, in particular the EOS 
and mass ratio, may lead to quantitative differences. The 
inclusion of microphysical aspects, neutrinos and mag¬ 
netic fields may change the picture [H], but because 
the mechanisms producing mass ejection described here 
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0 1000 2000 3000 4000 5000 

time I/Mq 

FIG. 15: Rest-mass conservation for MSlb-100150 and res¬ 
olution study. The plot shows results using resolutions R1 
(blue), R2 (red), R3 (black), and runs with and without the 
C step. Top: rest-mass; Middle: error of the rest-mass con¬ 
servation; Bottom: maximum density Pmax{t) normalized by 
the initial maximum density pmax{t = 0) 

operate on very short timescales of a few milliseconds 
during the merger, we expect differences only on longer 
timescales. 


IX. CONCLUSION 

In this work, we have investigated the merger remnant 
of neutron star binaries using ab initio numerical relativ¬ 
ity simulations which employ a conservative algorithm 
for the adaptive mesh refinement (AMR) technique. Our 
results are summarized in the following. 

(i) We have presented a new implementation of the 
Berger-Collela mesh refinement algorithm in the BAM 
code. The algorithm has been extensively tested in sin¬ 
gle star spacetimes focusing on its performances when 
combined with different reconstruction and prolongation 
operators and a standard artificial atmosphere treatment 
for the vacuum regions. 

The use of a correction step in the AMR algorithm 
significantly improves rest-mass conservation. In all our 
tests we found an improvement of at least a factor of 
^ 10 up to a factor of ~ 10®. However, mass conserva¬ 
tion depends also on the atmosphere parameters. Typi¬ 
cally, smaller atmosphere levels led to smaller violations. 


The choice of the restriction/prolongation operators can 
be delicate as well. The best mass conservation was ob¬ 
tained, for most of the cases, using the average restriction 
and a 2nd-order ENO prolongation. 

(ii) We have applied the conservative AMR in neutron 
star mergers simulations and focused on the study of the 
merger remnant. We considered initial binary configura¬ 
tions with different EOS, binary mass 2.7 Mq, and two 
mass ratios g = 1,1.16. Very similar simulations where 
performed in [55] . We studied the dependence of the 
merger outcome as a function of the EOS and q. For 
M = 2.7 Mq a massive differentially rotating object is 
produced, the properties of which mostly dependent on 
the EOS. Stiffer EOSs produce more stable remnants, 
and eventually stable objects (MNS) in cases the total 
rest-mass is less than the one supported by a spherical 
configuration with the same EOS. Softer EOSs produce a 
hyper-massive neutron star (HMNS) which collapses on 
dynamical timescales. The HMNS collapses to a black 
hole with mass Mbh 2.4 — 2.5 Mq and dimension¬ 
less spin ^ 0.58 — 0.64. An accretion disk of rest-mass 
Aldisk ~ 0.05 — 0.2 Mq and a radius of ^ 40km is ob¬ 
served. 

All the simulations were computed with and without 
conservative AMR. The conservative algorithm typically 
improved rest-mass conservation by a factor of ~ 5, de¬ 
pending on the specific resolution and binary configura¬ 
tion. At the resolutions employed, rest-mass violations 
can lead to inaccuracies in the collapse time, and to sys¬ 
tematic errors regarding the mass of the accretion disk, 
and the black hole mass/spin. Differences in the ejecta 
are also observed, although no general trend could be 
identified. Our results indicate that the use of conserva¬ 
tive AMR is desirable and recommended in postmerger 
simulations. 

We studied dynamical mass ejection and found that 
a total rest-mass of about Mejecta 10“® — 10“^ Mq 
becomes unbound during merger with kinetic energy 
Tejecta ~ 10“^ ~ 10®°erg. The amount of ejected ma¬ 
terial depends on the EOS and on the mass ratio q. For 
q = 1 larger ejecta are observed for softer EOSs. For 
a given EOS, larger q gives larger ejecta. Overall, our 
results agree with those of [SU |53| . Even with conser¬ 
vative AMR the computation of ejecta is challenging for 
numerical relativity grid-based codes, at least when the 
moving-box algorithm with nested boxes centered on the 
stars is used. Conceivably, a local AMR strategy that 
tracks the ejecta could be advantageous. In order to ob¬ 
tain the best performance one needs to carefully set and 
experiment with the atmosphere parameters. 

(iii) As a new application we have performed, for the 
first time, a simulation of a g = 1.5 configuration with the 
stiff EOS MS lb. Mass-ratio g = 1.5 is the largest mass 
ratio simulated so far in numerical relativity, simulated 
in for a very soft EOS. Here, we considered the very 
stiff EOS MSlb; the two stars have masses I.OOMq and 
I. 5 OM 0 . 

During merger the companion (less massive star) is 
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FIG. 16: Snapshots of the MSlb-100150-Rlc evolution on the {x,y,z = 0) and {x,y — 0,a;) planes for t = 2016Mq, 2794Mq 
(upper panels) and t = 2961M0, SOSIMq (bottom panels). The density p is plotted in logarithmic scale with white con¬ 
tours shown cit p — (10“^, 10“®, 10“®, 10“'*, 10“®), the ejecta are colored red (or black for better readability) at p = 
(IQ-io^ 10“®, 10“®, 10“^, 10“®), the velocity v* is visualized by black arrows. The logarithm of the entropy indicator logj^Q S is 
presented according to the color bar. 
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FIG. 17: Mass, average velocities and geometry of the ejected 
mass MSlb-100150. Top: Ejecta mass for all resolutions. 
Middle: (u)p.2 on the (a;, y) plane and the {x, z) plane. Bot¬ 
tom: (p), (z). The middle and bottom panel are restricted to 
highest and lowest resolution only for better readability. 



0 1000 2000 3000 4000 5000 6000 


time u/Mq 

FIG. 18: SLy-125145-R2 configuration for different grid se¬ 
tups. Without the G step, a change in resolution and box 
size of 3% has a large influence on the lifetime of the HMNS. 


electromagnetic (and neutrino) signals. They should be 
investigated in the future in more detail including mag¬ 
netic helds, microphysics and radiation transport in the 
simulations. 


strongly deformed by the tidal held of the primary and 
develops a tidal tail. The hnal merger remnant is com¬ 
posed by a high density hot core surrounded by a thick 
accretion disk of rest-mass ~ 0.3 Mq and of radius 
~ 35Mq ~ 55km (see Fig. 141. The remnant is not 
expected to collapse since the total binary rest-mass is 
smaller than the maximum rest-mass supported by a 
spherical conhguration. 


The MSlb-100150 conhguration has the largest 
amount of ejected rest-mass in our sample, Mejecta 
0.03 Mq . Ejecta mainly originate from the tidal tail; den¬ 
sity layers of order p ^ 10“® — 10“^ ^ 10® — 10^° g cm“® 
are accelerated up to ^ 0.3 and become unbound. Most 
of the unbound mass is ejected in a time window of a 
few milliseconds around the moment of merger, tmrg- 
We identihed two mechanisms for the ejecta emission: 
(i) the torque exerted by the central two-cores structure 
on the tidal tail; and (ii) shocks waves generated between 
the two MNS cores. The geometry of the emission is 
anisotropic. Although conhguration details may lead to 
some quantitative differences, we suggest that our qual¬ 
itative picture is rather robust and captures accurately 
the short timescale dynamics of the ejecta. 

We believe conhgurations like MSlb-100150 are astro- 
physically plausible and potentially relevant for strong 
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Appendix A: Box-sizes in BNS simulations 


In this appendix we investigate the influence of the box 
settings on the HMNS dynamics. In |53] we have exper¬ 
imented with the box sizes in an attempt of improving 
rest-mass conservation in the postmerger phase without 
using a conservative algorithm. We briefly compare here 
the two approaches. 

Focusing on SLy-125I45-R2, we consider runs with the 
different grid setting of Tab. IV In the grid setup R2[cn]l 
and R2[cn]2, the number of points per direction is kept 
fixed but the resolution is slightly changed in order to 
increase the box size. Changing the resolution has two 
competitive effects. On the one hand, the merger rem¬ 
nant is better resolved with a smaller grid spacing; but 
on the other hand the box size decreases and more matter 
can crosses refinement boundaries. 


Figure shows the central density and the gravita¬ 
tional wave signal. If no C step is employed, we observe 
a large shift 700Mq ^ 3.5 ms) in the collapse time. 
As an effect of the non-conservative AMR, the total mass 
increases and the system collapses earlier. In case the C 
step is applied a smaller shift of about ~ 300Mq ~ 1.5 ms 
is observed. This is possibly due to a similar effect as 
above, but of reduced magnitude and possibly due to the 
different resolution. No visible differences are observed 
in the GW signal instead. 

Increasing the box size while maintaining the same 
resolution increases the computational cost significantly, 
~ n^. On the other hand the computational overhead 
due to the C step amount to a maximum of 10% in sim¬ 
ulation speed, in the cases where the mask for the child 
and parent cells have to be computed often. We conclude 
the conservative AMR is a better approach. 
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